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We study high- field magnon dynamics and examine the dynamical structure factor in the quasi- 2D 
tetragonal Heisenberg antiferromagnet with inter layer coupling corresponding to realistic materials. 
Within spin- wave theory, we show that a non-zero interlayer coupling mitigates singular corrections 
to the excitation spectrum occurring in the high-field regime that would otherwise require a self- 
consistent approach beyond the 1/ S approximation. For the fields between the threshold for decays 
and saturation field we observe widening of the two-magnon sidebands with significant shifting of 
the spectral weight away from the quasiparticle peak. We find spectrum broadening throughout 
large regions of the Brillouin zone, dramatic redistributions of spectral weight to the two-magnon 
continuum, two-peak structures and other features clearly unlike conventional single-particle peaks. 
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I. INTRODUCTION 

Studies of quantum antiferromagnets are prominent in 
magnetism, elucidating the role of symmetry, fluctua- 
tions and dimensionality in the ground state of many- 
body systems, revealing complex dynamical properties 
of spin excitations and allowing detailed comparison be- 
tween theory and experiments. ^^^ A better understand- 
ing of many phenomena, from superfluidity to complex 
quantum phases may also be achieved through com- 
parative investigations of systems of strongly interact- 
ing bosons and frustrated antiferromagnets.^^ Recent 
developments in the synthesis of molecular based an- 
tiferromagnets with moderate exchange constanti^^^^ 
has opened previously unreachable high m agnetic field 
regime to experimental invest igat ions .^^^'^ Furthermore, 
the recent completion of neutron scattering instruments 
with enhanced resolution^^ ^^ provides an opportunity 
to study the dynamics of quantum antiferromagnets in a 
much wider momentum-energy space and under the in- 
fluence of applied magnetic field. 

In a collinear antiferromagnet, e.g., a square-lattice 
Heisenberg antiferromagnet, magnetic field induces a 
non-collinearity of the magnetic order, similar to the ef- 
fect of geometric frustration due to competing interac- 
tions, present in some other spin systems .^^^^ Gener- 
ally, non-collinearity leads to an enhanced interaction 
among spin excitations^^^ and results in stark differ- 
ences from conventional theory for magnons in the high- 
field regime .^^^^'^ In increasing field, spins gradually cant 
toward the field direction until they reach the satura- 
tion field, Hs^ where quantum fluctuations are fully sup- 
pressed and ferromagnetic alignment achieved. For fields 
above a certain threshold value, H* ^ but below satura- 
tion, coupling of the transverse and longitudinal spin 
fluctuations provides a channel for decays, through cu- 
bic terms in the spin-wave expansion. In this regime. 



magnons are predicted to be strongly damped, resulting 
in the loss of well-defined quasiparticle peaks.^^ 

Magnon decays and spectrum broadening above the 
threshold field in the purely two-dimensional square- 
lattice antiferromagnet have been confirmed via quan- 
tum Monte Carlo^^ (QMC) and exact diagonalizatioiP^ 
numerical studies. In addition, QMC has revealed a 
non-trivial redistribution of spectral weight resulting in 
non-Lorentzian "double peak" features in the dynamical 
structure factor, also previously observed in the origi- 
nal analytical study. ^^ The inelastic neutron-scattering 
experiment on the spin-5/2 material Ba2MnGe2 07 are 
also indicative of the field-induced magnon decays,!^ al- 
though not fully conclusive. ^^ It can be argued that re- 
cent thermal conductivity experiments in the so-called 
Bose-Einstein condensed magnets that observed suppres- 
sion of heat current in the vicinity of critical fields^^ are 
related to magnon decay dynamics. ^ Curr ent advances 
in neutron scattering instrumentation^^ ^^ open avenues 
to search for decays in large portions of the (k, uo) space. 
This is a particularly important issue in the case of spin- 
1/2 systems for which a comparison of experimental find- 
ings with existing theories is still missing. 

In this work, we extend the previous work of three of 
us and provide a theoretical investigation of high-field 
dynamics in the quasi-2D tetragonal 5 = 1/2 Heisenberg 
antiferromagnet with interlayer coupling corresponding 
to realistic materials. In particular, within spin-wave 
theory in the 1/S approximation, we obtain quantitative 
predictions for the dynamical structure factor, 6'(k,co'), 
the quantity measured in inelastic neutron scattering ex- 
periments, for a representative interlayer coupling ratio 
J7J = 0.2, relevant, e.g., for (SCAPjsCuCU,^^ along sev- 
eral representative paths in the Brillouin zone and for 
magnetic fields of H = 0.90Hs and H = 0.9bHs. 

Within the framework of spin-wave theory, 3D cou- 
pling also helps with the following technical issues. First, 



in the case of the 2D, 5* = 1/2 square-lattice antiferromag- 
net in high fields, calculations without the self-consistent 
treatment of cubic magnon interactions lead to renor- 
malized quasiparticle peaks that unrealistically escape 
the (unrenormalized) two-magnon continuum, preclud- 
ing an accurate determination of the dynamical struc- 
ture factor.^^ Second, for the same 2D case, the stan- 
dard 1/S expansion for the magnon spectrum breaks 
down in the high-field regime, owing to the transfer of 
the van Hove singularities from the two-magnon contin- 
uum to the o ne-magnon mode due to a coupling between 
the two.^^^ To regularize such unphysical singularities, 
two self-consistent schemes were developed in Refs. 21 
and [22l While the former approaclP^ (\i(\ take into ac- 
count the renormalization of the spectrum and, as a con- 
sequence, the spectral weight redistribution, it was only 
partially self-consistent. The latter, on the other hand, 
ignored the real part of the spectrum renormalization,^ 
essentially enforcing Lorentzian shapes of the quasipar- 
ticle peaks and excluding more complex profiles such as 
"double-peak" features . ^^^^ While this latter approach 
is suitable for spins S* > 1,^^ it cannot be deemed satis- 
factory in the case of 5 = 1/2. 



The aim of this work is to demonstrate that a non-zero 
3D interlayer coupling largely mitigates the singularities 
of high-field corrections and thus provides a physical rep- 
resentation of the excitation spectrum without the use of 
self-consistency. While the quasiparticle energy shift and 
corresponding escape from the "bare" continuum are still 
present to some degree within this approach, they are 
considerably weaker than in the purely 2D case, allowing 
for a detailed analysis of the effects of broadening and 
weight redistribution in the spectrum. Altogether, the 
relatively simple 1/S approximation requires consider- 
ably less computation than self-consistent methods, fol- 
lows the spirit of the regular spin-wave expansion and 
allows quantitative examination of the dynamical struc- 
ture factor. This approach is applicable to all spin values 
including the spin- 1/2 case, which is the focus of our 
study. Since all real antiferromagnets are at best quasi- 
two-dimensional, our analysis is relevant to most realistic 
materials, even those with weak interlayer coupling. Re- 
sults of the present work therefore consist in detailed and 
clear predictions for the dynamical structure factor to 
guide experiments and a comparison with self-consistent 
techniques for the purely 2D case. 



The paper is organized as follows: Section II contains 
an overview of the spin- wave formalism, with full de- 
tails in Appendix A. Following in Section III is a brief 
discussion of decay conditions and singularities. Section 
IV contains the dynamical structure factor calculated in 
the 1/5* approximation for non-zero interlayer coupling. 
Section V contains our conclusions. Appendix A con- 
tains details of the 1/S spin- wave theory for the quasi-2D 
tetragonal Heisenberg antiferromagnet in external field. 
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FIG. 1. (Color online) Self-energies in the lowest 1/S order 
obtained from decay (left) and source (right) interactions. 



II. INTERACTING SPIN WAVES IN FIELD 

We begin with the Heisenberg Hamiltonian of nearest- 
neighbor interacting spins on a tetragonal lattice in the 
presence of an applied magnetic field directed along zq 
axis in the laboratory reference frame. 
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where in-plane coupling is Jij = J and the interplane one 
is Jij = J' = a J and we assume < a < 1. With the details 
of the technical approach explicated in Appendix [A| and 
Refs. [2T| |2?, and '30, we summarize here the key steps of 
the spin-wave theory approach to this problem. 

First, we identify the canted spin configuration in the 
equivalent classical spin model. Then we quantize the 
spin components in the rotating frame that aligns the 
local spin quantization axis on each site in the direc- 
tion given by such a classical configuration. Applica- 
tion of the standard Holstein-Primakoff transformation 
bosonizes spin operators. After subsequent Fourier trans- 
formation, the Hamiltonian is diagonalized via the Bo- 
goliubov transformation, yielding the Hamiltonian that 
can be written a^^ 

^ = E ~^^^iK + ^ E ^1 (k, q) {hi_^^^h\h^ + h.c.) (2) 
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where the ordering vector Q = (7r,7r,7r) enters the mo- 
mentum conservation condition because of the staggered 
canting of spins.l^ In this expression, e\^ = 6:k + fek, 
where £\^ is the "bare" magnon dispersion given by lin- 
ear spin- wave theory and 5£\^ contains 1/S corrections 
from angle renormalization and Hartree-Fock decoupling 
of cubic and quartic perturbations, respectively. Ellipses 
stand for higher-order terms in the 1/S expansion that 
are neglected in our approximation. The three-boson 
terms are decay and source vertices that are responsible 
for the anomalous dynamics in the high-field regime. The 
dynamical properties of the system are obtained from the 
interacting magnon Green's function, defined as 



G-^(k,a;) 



■£k-Si(k,cj)-S2(k,cj) (3) 



where Si^2(k, cj) are the decay and source self-energies 
presented in Fig. [l] and obtained from the second-order 



treatment of Eq. ^ . Their explicit forms are 

E,(k,c.) = lV ^i(^ 

q 



a = 0.2 



^2{K0J) 



q — ^k-q+Q + ^0 

|^2(k,q)P 



2^^ uj- 
q 



^q + ^k+q-Q — ^0 



(4) 
(5) 



with expressions for vertices $i and ^2 given in Ap- 
pendix |A] In contrast to H = case where cj-dependent 
magnon interactions beyond Hartree-Fock have l/S"^ 
smahness, they readily occur in 1/S order for H y^ 0^ 
when the coupling between longitudinal and transverse 
modes renders three-boson vertices nonzero. Above the 
threshold field for magnon instability, the self-energy in 
Q also acquires an imaginary component, signifying the 
occurrence of spontaneous decays. 



III. KINEMATICS, SINGULARITIES AND 
INTERLAYER COUPLING 

A. Decay boundaries 

Although the single- and two-magnon continuum exci- 
tations are coupled directly by virtue of the cubic terms 
in Eq. ^ at any H > 0^ magnon decays only occur 
above a finite threshold field H^ . This is due to restric- 
tions provided by the kinematic conditions, i.e., energy 
and momentum conservation, that have to be satisfied in 
each elementary decay process 



^k — ^q + ^k-q+Q, 



(6) 



where the ordering vector in the momentum conservation 
is, again, due to the staggered canting of spins in the field. 

With the detailed classification of possible solutions of 
Eq. ([6| given previously in Refs. l2Ql and [22l we highlight 
here two relevant results. First, it can be shown that 
within the Born approximation, that is, neglecting energy 
renormalization of the spectrum Sk from the linear spin- 
wave theory result, the value of the threshold field H^ is 
independent of the value of the interlayer coupling a and 
is the same as in the square lattice case: H"^ ^0.76Hs. 

Second, for a given field H > H*, within the same 
Born approximation that neglects cascade decays, there 
exists a sharp boundary in the momentum space sepa- 
rating the region where magnons are stable, i.e., cannot 
satisfy Eq. (|6|, from the decay region where they are un- 
stable towards decays. An example of such decay bound- 
aries in a 3D setting is shown in Fig. |2] in one octant of 
the tetragonal-lattice Brillouin zone for a = 0.2 and sev- 
eral fields. ^^ One can see that the decay region expands 
with increasing field, covering large portions of the Bril- 
louin zone. At high fields, decay regions extend beyond 
the octant and overlap; at H ^ 0.9 Hg most magnons are 
unstable already in the Born approximation. 

Taking into account the finite lifetime of the de- 
cay products in the self-consistent treatment generally 




FIG. 2. (Color online) Magnon decay regions in one octant of 
the Brillouin zone for several fields. Decays are possible within 
the portion of the Brillouin zone containing the Mtt point. 
At high fields, decay regions extend beyond the octant and 
overlap, facilitating larger phase space for decays for magnons 
within the decay region. 



lead s to b lurring away the decay boundaries discussed 
j^g^g [21122] Nevertheless, it is still important to consider 
them not only because the decays within the Born de- 
cay region typically remain more intense as they oc- 
cur in the lower-order process, but also because such 
boundaries correspond to singularities in the spectrum 
discussed next.^ 



B. On-shell approximation 

Within the standard spin- wave expansion, the on-shell 
approximation for the magnon spectrum consists of set- 
ting uj = £\^ in the self-energies in Eqs. Q, ([s]), leading 
to the renormalized spectrum 



^k = ek + Kei;i(k,£k) + i;2(k,£:k), 



(7) 



and the decay rate 

Tk = ^ E l*i(k,q)l'5(<rk -e^- Sk-q+q). (8) 



We note that this approach implicitly suggests that the 
dynamical response can be well approximated by the 
quasiparticle peaks at renormalized energies £k with 
lorentzian broadening defined by Fk- As was demon- 
strated in Ref. [21] and [22] for the 2D square-lattice 
and in Ref. [20] for the triangular-lattice cases, this ap- 
proach, which is normally very accurate even quantita- 
tively, breaks down rather dramatically, showing diver- 
gences in £k or Fk for the momenta belonging to various 




FIG. 3. (Color online) Magnon energy and decay rate for 
a — 0.2 and for the square-lattice case [a — 0) for H — 0.9Hs 
along the selected path in the Brillouin zone (see Fig. ItI. 
Upper and lower panels show the renormalized spin- wave en- 
ergy £k and the decay rate Fk in the on-shell approximation, 
Eqs. ([7| and ([s}, respectively. Computation was performed 
with Monte-Carlo numerical integration with 5 • 10^ steps. 



contours in k- space. We reproduce some of them in our 
Fig. [3J shown by dashed lines. These singularities were 
understood as coming from the van Hove singularities in 
the two-magnon continuum through coupling with the 
single-magnon branch. One can see in the expression for 
Fk in Eq. ([8| that once the continuum and the single- 
magnon branch overlap in some portion of the Brillouin 
zone, ^k is effectively exploring the two-magnon density 
of states as a function of k. If a singularity in the latter 
is met, it will be reflected in Fk and, by the Kramers- 
Kronig relation, in £k- 

Generally, there are two types of singularities that can 
occur, associated with the minima (maxima) and saddle 
points of the two-magnon continuum, respectively. The 
decay boundaries discussed previously are necessarily the 
surfaces of such singularities because they must corre- 
spond to the intersection with a minimum of the two- 
magnon continuum, forming the locus of points where 
decay conditions are first met. These decay threshold sin- 
gularities correspond to the step-like behavior in Fk and 
the logarithmic divergence in £k-^^ ^^ For the k-points 
already inside the decay region, there is a possibility of 
meeting a saddle point of the continuum, in which case 
the real part of the spectrum Sk experiences a jump and 
Fk a logarithmic singularity, see the 2D data in Fig. [3] 

Essential for the present study is the change in the 
behavior of these singularities for a moderate interlayer 
coupling. One can expect from the density of states 
argument that both types of singularities must become 
weaker and change from logarithmic/step-like to square- 
root-like ,1^ resulting in softening of the singular behavior 
of Fk and £k- This effect is demonstrated in Fig. [3] for 
several directions in the 3D Brillouin zone compared with 
equivalent cuts for the square-lattice case. The full ac- 
count of analytical and numerical details of the on-shell 



calculations is given in Appendix A. The definitions con- 
cerning the Brillouin zone path in Fig. [3] can be found 
in Fig. 7. It is rather remarkable that the singularity 
of the saddle-point type is completely wiped out by the 
moderate interplane coupling, see panel F — X in Fig. [3] 
While the decay threshold singularities associated with 
the boundaries in Fig. [2J are not fully alleviated, they 
are considerably diminished. For both types of singu- 
larities, further increase of the interplane coupling leads 
to gradual changes in the on-shell spectrum and decay 
rates; singularities shown in Fig. |3] for a = 0.2 remain 
largely unchanged. 

While singularities are diminished, the magnon damp- 
ing remains considerable. Since both the decay bound- 
aries and Fk are obtained within the Born approxima- 
tion, the decay rate in Fig. |3] is non-zero within the cor- 
responding decay region outlined in Fig. [2] for H = 0.9Hs, 
demonstrating that magnons are unstable in most of 
the Brillouin zone, while the Goldstone mode at M^^ 
[k = (tt, TT, tt)] and the uniform precession mode at F 
[k= (0, 0, 0)] remain well-defined. 

Previous works on regularization of singularities i n the 
spin-wave spectrum in the 2D square-lattice case^^^ 
were based on partial dressings of magnon Green's func- 
tions in the one-loop diagrams of Fig. [l] Generally, the 
difficulty of a consistent implementation of such schemes 
is the opening of an unphysical gap in the acoustic 
branch. The latter problem was avoided in Ref. |22] by 
performing self-consistency only in the imaginary part 
of the magnon self-energy, applicable for larger S > 1 
where the real part of renormalization can be deemed 
small. The downside of this approach is that it remains 
essentially on-shell, enforcing Lorentzian shapes of the 
quasiparticle peaks. In a sense, what is being argued by 
our Fig. [3] is that introducing finite interplanar coupling 
represents an alternative to the self-consistency schemes 
in regularizing singularities, without restrictions on the 
spectral shapes and at a fraction of the computational 
cost. 



C. Spectral function 

The one-loop self-energies 111^2 (k,cj) in Eqs. Q and 
([5]) originate from the coupling to the two-magnon con- 
tinuum. Because of their cj-dependence, one can obtain 
significantly richer dynamical information, beyond the 
quasiparticle pole-like state expected from the on-shell 
approach above, by considering the diagonal component 
of the spectral function. 



A(k,cj) 



1 



ImG(k,cj), 



(9) 



where G(k, cj) is the interacting magnon Green's func- 
tion defined in Eq. (|3|. Because of the interactions, the 
spectral function is expected to exhibit an incoherent 
component which should reflect the two-magnon contin- 
uum states in addition to a quasiparticle peak that may 



be broadened. The spectral function is directly related 
to the dynamical structure factor S'(k,a;) measured in 
neutron-scattering experiments, which we discuss in de- 
tail in the next Section. Strictly speaking, the considera- 
tion of A(k, u) goes beyond the 1/S expansion as a; 7^ £k, 
yet it is free from the complications due to higher-order 
diagrams mentioned above. 

We should note that in the pure 2D, S = 1/2 case, non- 
selfconsistent calculation of A(k,uj) does exhibit incoher- 
ent subbands,^but suffers from the high density of states 
of the two-magnon continuum associated with the thresh- 
old singularities considered above. Because of level re- 
pulsion, the renormalized quasiparticle peaks escape the 
unrenormalized two-magnon continuum and preclude an 
accurate determination of the dynamical structure factor, 
except for fields in close vicinity of the saturation field 
where the renormalization becomes weaker. However, 
this is not the case in the triangular-lattice antiferromag- 
net in zero field, ^^ where the thresholds are controlled by 
the emission of the Goldstone magnons associated with 
much smaller density of states. 

Similar to the weakening of the threshold singulari- 
ties in 3D, one can expect weaker repulsion between the 
single- magnon branch and the minima of the two-magnon 
continuum. As will be shown in the next Section, while 
the quasiparticle energy shift and corresponding escape 
from the continuum are still present to some degree, they 
are considerably smaller than in the 2D case. Therefore, 
once again, introducing finite interplanar coupling pro- 
vides an alternative to the self-consistency schemes, al- 
lowing for a detailed analysis of the effects of broadening 
and weight redistribution in the spectrum. 



IV. DYNAMICAL STRUCTURE FACTOR 

The inelastic neutron-scattering cross section is pro- 
portional to a linear combination of the components of 
the spin-spin dynamical correlation functioiP^ 



PC 
J — C 



dt 

2^ 



S^{t)S\)e 



Aojt 



(10) 



where a, /3 = {xq, ^o? ^0} span the Cartesian directions of 
the laboratory frame. For simplicity, in the following dis- 
cussion we ignore additional momentum-dependent po- 
larization factors from experimental details, on which co- 
efficients of such linear combination can depend on, and 
show our results for what we will refer to as the "full" 
dynamical structure factor 

^(k, u) = 5^°^° (k, uj) + sy^y^ (k, uj) + s'^'^ (k, cj), (11) 

as well as its component perpendicular to the external 
field, the latter is directed along the zo-axis as before, 

^^(k,cj) = ^^°^o(k,cj) + ^^°^°(k,cj). (12) 

The choice of these particular examples is both illustra- 
tive and general, as it exposes all the essential compo- 
nents of 6'(k, cj). 



The dynamical structure factor is naturally written in 
a laboratory reference frame while the magnon opera- 
tors were introduced in the rotating frame that aligns 
local spin quantization axis on each site in the direction 
given by a classical configuration of spins canted in a 
field. Thus, the relation of the magnon spectral func- 
tion, A(k, cj) of Eq. ([9|, to 5'(k, cj) for canted spin struc- 
tures is via a two-step transformation: rotation of spins 
from laboratory to local frame and bosonization of spin 
operators.-^ In the lowest l/^-order, this procedure be- 
comes rather straightforward since most complications 
such as off-diagonal terms in the spin Green's functions 
can be justifiably dropped as they only contribute in the 
order 1/S'^ or higher. With that, performing the second 
transformation first, components of the dynamical struc- 
ture factor in the local quantization frame are readily 
given by^ 

5^^(k,cj) ^ 7r5'A+ {u^, + v^,fA{\i,uJ), 

Syy{\i,uj) ^ ttSA- (iik - vi,fA{k,uj), (13) 



where x^y^z now refer to local axes, u\^ and v\^ are pa- 
rameters of the Bogoliubov transformation, and K± = 
(1 — (2n^5)/2S) contain the 1/S corrections from the 
Hartree-Fock spin reduction factors, see Appendix A. 

Strictly speaking, in these expressions for the diago- 
nal terms of S^^ we have kept contributions in excess of 
the leading order of 1/S expansion, as both the Hartree- 
Fock corrections to S^^ and S^^ as well as the S^^ term 
itself are of higher order. In addition to being easy to 
include them into our consideration, they also offer us 
an opportunity to demonstrate their relative unimpor- 
tance. Thus, the longitudinal component {S^^) in the 
given order contains a "direct" contribution of the two- 
magnon continuum to the structure factor. Apart from 
being broadly distributed over the (k, uj) space compared 
to a more sharply concentrated 74(k, cj), its contribution 
to the scattering turned out to be diminishingly small 
compared to the two transverse components S^^^^^^ in 
the considered high- field regime. Besides being an effect 
of the higher 1/S order, this is also due to its dependence 
on v'^ (cxcos^ 6><C 1), where 6 is the spins' canting angle 
which tends to 7r/2 as the saturation field is approached. 

Completing the connection of the local form of cor- 
relation functions to the laboratory ones gives the final 
relation of A(k, cj) to 6'(k,cj) 

5^°^°(k,cj) ^sin26>5^^(k,cj)+cos26>5^^(k-Q,a;), 
^^o^°(k,a;)^cos26>5^^(k-Q,a;)+sin26>5^^(k,a;), 

sy^y^{k^uj) = syy{k,uj) , (i4) 

where the cross-terms in spin-spin correlation function 
{S^^ and S^^) in the r.h.s. were omitted under the same 
premise of not contributing in the lowest order. These 
last expressions demonstrate that due to external field, 
spin-canting redistributes the spectral weight over two 
transverse modes, 1^^^ referred to as the "in-plane" and 




0.5 1-0 

u;/(2 + a) 



2.0 




t^ 



0.0 



0.5 1-0 

u;/(2 + a) 



FIG. 4. (Color online) Profiles of the perpendicular component of the dynamical structure factor S±(k,(jj), Eq. (12), vs cj for 
selected k-points along the main diagonal of the Brillouin zone (from F to M^r), panels (a) and (b), and along the diagonal 
in the kz — ti plane (from F^^ to M^r), (c) and (d), for fields 0.9Hs and 0.95ifs, and S — 1/2. For high-symmetry k-points 
notations, see Figs. [2] and [7| Dashed lines indicates the linear spin- wave dispersion £k along the same paths. 



the "out-of-plane" modes, corresponding to the momen- 
tum k and to the momentum k — Q (and to fluctuations 
in and out of the plane perpendicular to applied field), re- 
spectively. As can be expected form the pre-factor cos^ 
in Eq. ([l4|, strong magnetic field suppresses the out- 
of-plane mode until it disappear entirely at Hg. In the 
5'(k,co') plots that will follow, the out-of-plane contribu- 
tion appears as a "shadow" of the main signal, shifted by 
the ordering vector. 



From Eqs. (14) and (13) one can see that our choice 



of the perpendicular component of the structure factor 
S'^(k, uj) in Eq. ( 12 ) as a representative example becomes 
particularly simple. Since it is not contaminated by the 
k — Q shadow component and the contribution of the 
longitudinal 5'^^(k, cj) to it is exceedingly small in the 
considered field range, it is closely approximated by 



5i(k,w)«7r5/kA(k,w), 



(15) 



where /k = [sin^^A+ (iik+i'k)^+A_ (i^k — '^k)^] is the k- 
and field-dependent intensity factor. 

Our Fig. |4] presents the perpendicular component of 
the dynamical structure factor 6'^(k, uo) for 6* = 1/2, fields 
of 0.9i^s and 0.95i^s, and for several k-points along two 
representative directions: main diagonal of the Brillouin 
zone (from F to M^^), Fig. [4]^a) and (b), and diagonal 
of the kz = TT plane (from F^^ to Mt^), Fig. ^c) and (d), 
respectively. To obtain these profiles, integrals in the self- 
energies in Eqs. (p| and ([5| and in S'^^(k, cj) in Eq. (Il3| 



were computed using Mathematica via adaptive quasi- 
Monte- Carlo without symbolic preprocessing, over 5-10^ 
points, maximum recursion of 10^, and an accuracy goal 
of 4 digits. Each presented k-cut contains 400 points in uj 
and an artificial broadening d/{2-\-a) = 10~^ J was used. 
As we discussed, contribution of the longitudinal term 
S^^ is not visible on the scale of the plots. Dashed lines 
show linear spin- wave dispersions for the given paths. 

These 6'x(k, u) profiles demonstrate that the spectrum 
broadening due to magnon decays in the high- field regime 
is found abundantly throughout the Brillouin zone. In 
addition, clear double-peak structures, similar to that 
seen in the self-consistent spin-wave calculation^^ and 
emphasized in the QMC stud}'^^ for the 2D square-lattice 
case, occur along the main diagonal. Fig. Qa), (b). At 
H = O.dOHg^ some sharp magnon peaks remain in the 
vicinity of the F-point, in agreement with the decay 
boundaries in Fig. |2] Peaks also get narrower at the 
approach of the Goldstone M^^ point. Although there 
seems to be some sharpness in the structure of S±{'k,uj) 
data for H = 0.90Hs in the middle of the diagonals, these 
are not associated with quasiparticle peaks, but with the 
remaining singular behavior of the one-loop self-energies. 

A noteworthy feature of the shown results is a strong 
field-evolution of 6'^(k, a;), demonstrating a rather dra- 
matic redistribution of spectral weight between the 
broadened single-particle and incoherent part of the spec- 
trum, related to the two-magnon continuum. Another 




FIG. 5. (Color online) Intensity plots of the dynamical structure factor 5'(k, cj), Eq. (11), for fields 0.9Hs and 0.95ifs, and 
S = 1/2, for the k-path in the /c^ =0 plane, see Fig. JTl Gaussian convolution was performed in u with a cr-width of ^cj^O.lJ. 



important observation is the relative smallness of the en- 
ergy shift of the spectral- function leading edge in uj com- 
pared to the linear spin- wave energy Sk, shown by the 
dashed lines. This, together with a close similarity of 
the shown non-selfconsistent results to the self-consistent 
ones in the 2D square-lattice case,^ is an indication of 
the reliability of current approach, thus supporting our 
expectation that presented results should serve as a fair 
representation of the realistic structure factor. 

In Figs. [5] and [6} intensity plots of the full dynami- 
cal structure factor, 6'(k, cj) in Eq. (11), for S = 1/2, 
fields of 0.9Hs and 0.95J^s, and for several representa- 
tive paths in the Brillouin zone are shown. In Fig. |5l 
the k-path is entirely in the kz = plane and in Fig. |6] 
the momentum traces an inter-layer path, both shown 
in Fig. [7| As we have discussed above, in addition to 
the main features demonstrated in Fig. [4] for 6'^(k, cj) 
that are directly related to the magnon spectral function 
74(k, cj), the out-of-plane transverse mode contribution is 
visible as a shadow, shifted by the ordering vector. The 
longitudinal (S^^) contributions are hardly noticeable on 
the scale of the plots as before. In Fig. [6) proximity 
to the Goldstone mode at M^ causes unphysical diver- 
gences, amplified by the divergent intensity factor /k in 
Eq. (15), that necessitate removal of this region in the 
plot. 

Integrals in the self-energies for these plots were com- 
puted by the same method as in Fig. [4j over a maxium 
number of points of 2 • 10^, with the same accuracy goal 
of 4 digits. There are 150 points along each segment in 
the k-path, each with 350 points in u. All integrations 
were performed in parallel over 8 logical threads. The 
same upper cut-off in intensity has been used in both 
plots with additional red color fill up to the extent of the 
highest peak, whose value is shown in the color sidebars. 
Additionally, we have performed Gaussian convolution in 
the cj-direction with a a- width of Suj = 0.1J. This step 



is intended to mimic the effect of a realistic experimen- 
tal resolution, from which we can also draw quantitative 
predictions of the relative strength of the quasiparticle 
and incoherent part of the spectrum. 

Despite the provided broadening by the finite en- 
ergy resolution, complex spectral lineshapes, very much 
distinct from the conventional quasiparticle peaks, are 
clearly visible in Figs. |5] and |6] This demonstrates that 
the effects of spectral weight redistribution and broaden- 
ing due to spontaneous decays are substantial and should 
be readily observed in experiment. 

Since F, M, X, and X' points in Fig.[5]and F, X, and X^ 
points in Fig. |6] are outside of the Born approximation de- 
cay regions according to Fig. |2J spectrum in their vicini- 
ties exhibits well-defined quasiparticle peaks, broadened 
by our "instrumental" resolution. However, away from 
them, the structure factor demonstrates a variety of un- 
usual features, including the already discussed double- 
peak lineshape for FM direction in Fig. [5] and FMt^ di- 
rection in Fig. [6J the latter path also shown in S^Ck^uj) 
previously. 

Along some of the k-directions for H = 0.9 Hs in both 
Fig. [5] and [6) the renormalized quasiparticle peaks escape 
the unrenormalized two-magnon continuum and survive 
despite being formally within the decay region. Yet, these 
peaks are accompanied by continuum-like subbands that 
accumulate significant weight. Upon increase of the field, 
the continuum overtakes the single-particle branches in 
the large regions of the Brillouin zone, washing them out 
and creating more double-peak structures, e.g., very dis- 
tinctly along the XT (Fig. [5| and X^F (Fig. ^ direc- 
tions. There is a particularly spectacular advance of the 
continuum along the XX' line in Fig. [5J where magnon 
broadening is most intense with dramatic spectral weight 
transfer to high-energies. 

We note that the regions where decays are present in 
the on-shell solutions discussed in Sec. Ill often corre- 




FIG. 6. (Color online) Same as in Fig. [s] for a different path, see Fig. |7| The region near M^r in k and k — Q modes are cut 
due to spurious divergences, see text. Note the higher cutoff in intensity due to a formfactor divergent at Mtt -point. 



spond closely to the "washout" regions of the intensity 
plots, although the latter draw a much richer picture by 
exposing a complex cj-structure. The remaining 3D sin- 
gularities in the real part of the on-shell energy, exempli- 
fied in Fig. [3J indicate where the two-magnon continuum 
overlaps with the single-magnon branch, also in a quali- 
tative agreement with the intensity plots, thus providing 
a complimentary information. Previous examination of 
the 2D square-lattice case in Ref. 22 showed broaden- 
ing patterns similar to the ones in Figs. |5] and Fig. [6] for 
corresponding paths of the Brillouin zone, although the 
cj-structure differs significantly as the self-consistent ap- 
proach of that work has enforced the Lorentzian shapes 
of the spectral function. 




FIG. 7. (Color online) Brillouin zone octant showing the k- 
paths traversed in Figs. [5] and [6] Path in Fig. [5] lies in the 
fe = plane and goes along the dotted line: F ^ M ^ X ^ 
X' -^ F. Path in Fig. [6] goes between the kz = and kz = tt 
planes along the zig-zag solid line: X^^r ^ F ^ M^r -^ X. 



The one-loop approximation used in this work takes 
into account the real part of the self-energy corrections, 
allowing the renormalized single-magnon branch to es- 
cape some of the decay region. Therefore, the obtained 
dynamical structure factor is likely to overestimate the 
field values needed to reach a significant overlap between 
the two-magnon continuum and the single magnon mode 
throughout the Brillouin zone. 

Altogether, examination of the dynamical structure 
factor reveals rich features consisting of dramatic redistri- 
butions of magnon spectral weight. This ranges from the 
appearance of double-peak lineshapes to a full suppres- 
sion of quasiparticle peaks throughout large portions of 
the Brillouin zone. These features fully survive convolu- 
tion with moderate energy resolution and should there- 
fore be accessible to state-of-the-art neutron scattering 
experiments. 



V. CONCLUSIONS 

In this work, we have provided a theoretical consider- 
ation of high-field dynamics in the quasi-2D tetragonal 
S = l/2 Heisenberg antiferromagnet with interlayer cou- 
pling corresponding to realistic materials. We have pre- 
sented quantitative predictions for the dynamical struc- 
ture factor for a representative interlayer coupling ratio 
J' /J = 0.2 along several paths in the Brillouin zone. We 
have demonstrated within the framework of spin-wave 
theory that the finite 3D interlayer coupling, present to 
some degree in all realistic materials, largely mitigates 
singularities of high-field corrections that are complicat- 
ing similar analysis in the purely 2D case and necessitate 
a self-consistent regularization in the latter. 

We argue that introducing finite interplanar coupling 
effectively provides an alternative to the self-consistency 
schemes, allowing for a detailed analysis of the effects 



of broadening and weight redistribution in the spectrum 
without the use of self-consistency. This relatively simple 
procedure requires considerably less computation than 
self-consistent calculations, follows the spirit of the 1/ S 
spin-wave expansion and provides a physical represen- 
tation of the excitation spectrum. Our results there- 
fore consist in a detailed picture of the latter and offer 
clear predictions to guide experiments. Such a presenta- 
tion is increasingly valuable as state-of-the-art neutron- 
scattering instrumentation should allow the spin dynam- 
ics of real materials to be searched for the discussed 
unusual features with sufficient energy resolution and 
momentum-space coverage. 

The present analysis does not include the case of frus- 
trated interplanar exchanges for which enhanced quan- 
tum fluctuations^^ should result in decays of comparable 
magnitude to the purely 2D case. Additionally, ferromag- 
netic interlayer coupling, while experimentally appealing 
due to the reduction of the threshold field H*^^^ goes be- 
yond the scope of the present work. Nonetheless it is ex- 
pected that the inclusion of any type of increased dimen- 
sionality of sufficient strength will act to soften Van Hove 
singularities without inherently removing decay phenom- 
ena. 

The landscape of the high-field dynamical structure 
factor shown in this work is diverse and intriguing. Our 
results, consistent with prior numerical and analytical 
studies, provide a comprehensive illustration of its de- 
tails. Inelastic neutron-scattering measurements on suit- 
able quantum antiferromagnets will be important to elu- 
cidate the accuracy of the theoretical results presented 
in this work. 
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Appendix A: Details of Spin- Wave Theory 

The derivation of 1/S corrections to spin- wave the- 
ory in an applied magnetic field with non-zero interlayer 
antiferromagnetic exchange begins with the Heisenberg 
Hamiltonian of nearest-neighbor exchange and Zeeman 
contributions. 



n 



E 

(ij) 



JijSi-Sj H y ^ J j ° , 



(Al) 



Where the factor gjHB^ with g being the gyromagnetic 
ratio and jj^b the Bohr magneton, is included in the 



applied magnetic field H. We consider the Heisenberg 
antiferromagnet on the simple tetragonal lattice, with 
nearest-neighbor exchange in-plane coupling J and inter- 
layer coupling J', parametrized by a = jy J, < a < 1. 
The long-range antiferromagnetic order of spins canted 
by external field is described by the same ordering wave- 
vector Q = (tt, TT, tt) as in the collinear Neel order at 
H = 0. 

Spin components in the laboratory frame (xo,^o,^o) 
are related to that in rotating frame (x, y^ z) via the 
transformation 



^f = S^ sin - e^^-^^^f cos 6>, 

gxo = e^'^-'^^S^ COS O^Sf sine, 
qyo _ qy 



(A2) 



with 9 the canting angle from the spin-flop plane. The 
transformed Hamiltonian reads 



-e'^-'"' sin 20{SfS] - S^S^) (A3) 



iij) 



-HJ2 i^f sin - e**^'"' Sf cos 6') , 

i 

where (ij) designate a single counting of nearest neigh- 
bor exchanges. Using the usual Holstein-Primakoff 
transformation^^ and keeping up to quartic terms we ar- 
rive at the following terms of order (9(6'^~^/^): 



eo = Uq/N = - J6'M (2 + a) cos 2(9 + sin 6 



H 

'Js 



Hi = cos6> ^ e^^-^^ S^{H - AJSsinO (2 + a)) 

i 

I-L2 = 2^ iHsinO a\ai 

y^ Jij ( cos 20 (a\ai + o}-aj) 



S_ 
2" 



1^ {a\aj + a/jai) 



{a\a}- + ajGi 



n. 



E< 



,iQr, 



iJcos6' J(2 + a) sin 26* 

~^ 2 



(A4) 
(A5) 

(A6) 

(A7) 



:{a\a\ai -^ a\aiai) + sin 26> ^ J^^ (aj -^ ai)o}-aj 



H4 - 2^ Jv 



.j{i) 



■ cos^ 6 {[rii + rij) aiaj + h.c.) 



(A8) 



- sin^ (a\{ni + rij) aj -\- h.c.) — cos 2^ n^nj 



where j{i) refers to a site j that is nearest neighbor of i. 
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1. Harmonic approximation 

First, we minimize Eq to obtain the classical value of 
the canting angle 

H H 



smt 



4JS{2 + a) Hs 



(A9) 



where Hs =475(2 + a) is the saturation field where the 
spin structure reaches fully polarized state. Then, the 
Fourier transformation of I-L2 leads to the usual 

^^2 = XI (^k4^k - 2^k (4«-k + «ka-k) j , (AlO) 



where the coefficients A\^ and B\^ read 

Ak = 2JS(2 + a) (1 + sin^ 7k) 
5k = 2JS{2 + a) cos^ i9 7k 



(All) 



and where the lattice harmonics are defined as 7k = 
{ixy + ^lz)/{'^ + ^) with ^xy = (cos/ca; + COS ky) and 
7^ = cos kz . Upon the Bogolioubov transformation and 
substitution of the classical canting angle, we obtain the 
dispersion relation. 



£:k = 2J5(2 + a)cJk, 



(A12) 



cjk = \/(l + 7k)(l-cos2l9 7k) 
and the parameters of the transformation 



2 2 _ ^k i ^k 



'^^k'^k = - — . (A13) 
2sk 



2. Mean-Field decoupling 



default global adaptive integration methods via Mathe- 
matica with an accuracy goal of greater than 5 digits. 

Using the method of Oguchi^^ to treat the quartic 
term, we obtain the quadratic Hamiltonian, 



k 

where 



( ^k^-k + ^ktt-k j , (A15) 



SA 



[^^ = 2 J(2 + a 

, 7k 



A cos^ — n cos 20 — m sin^ 
(J cos^ — 2n sin^ 0) — jrn cos 20 



?(4) 



5Bi;^ = 2 J(2 + a) 

, 7k 



(A sin^ ( 



'0) 



(A16) 



{S sin^ -2n cos^ 0) + 7a cos 20 



using the definitions 



2m^ 



- arriz 



2 + a 



^^2A., + aA.^ (A17) 



7m = ^ ^ "^ , ^, 7A 



2 + a 



2+a ' ' 2+a 

The corresponding correction to the dispersion reads 
,(1) _ 9 7^9 ^ ^A (A18) 



fej^'^ = 2 J(2 + a) 

X — [(1 + 7k sin^ 0)8A^^^ - 7k cos^ 08B^ 

CJk 1- 



where A\^ a nd ^ k are dimensionless expressions in the 
brackets of (A16). 



The quartic Hamiltonian H4 contains four-boson terms 
that can be treated using mean-field Hartree-Fock decou- 
pling. To do so, we introduce the following averages. 



n= {a\ai) = ^vl , 

k 

3 = (aitti) = ^^Xk^k 



(A14) 



_ / t \ _ Y^ ^^y 2 

"mxy — [(^i^lj/xy — 2_^ ^~'^k 5 
k 

ruz = {alaj)z = '^lzvl , 

k 

^xy = {(^iCij)xy = ^ ^'^^k'^k 
k 

A^ = {aiaj)z = ^72^k^k • 



In zero magnetic field, anomalous averages vanish so that 
^ = '^xy = rriz = while the average n determines the 
sublattice magnetization associated with the linear spin- 
wave theory (S^) = S — n. Evaluation of the Hartree- 
Fock averages in finite magnetic field is performed using 



3. Angle Renormalization 

In contrast to the zero-magnetic field case, the 1/S ex- 
pansion of the present Hamiltonian contains cubic terms. 
The Hartree-Fock decoupling of the cubic term Hs is re- 
sponsible for a renormalization of the canting angle and 
yields 



(%) = V25J(2 + a) (A19) 

X sin 20{n — A — m) [ciq -\- clq) • 

The value of the renormalized angle is obtained from the 
cancellation of the linear term Hi + (H3) so that 



1 



n — m — A 



(A20) 



and the corresponding correction to the dispersion is 
therefore 



Se- 



(2) 



4 J(2 + a) sin^ 

(A-\-m — n) 
X- 

CJk 



(A21) 



1 - 7k - 7k cos^ / 
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Cubic Vertices 



Interaction vertices are obtained upon Fourier trans- 
formation of the fluctuating part of Hs , 



n3 = V2SJ{2^a) 

X sin 20 ^ 7k (^aj^a^ak-q+Q 

k,q 



") 



(A22) 



leading to the following cubic decay and source ver- 
tices, <l>3i(k,q) and <l>32(k,q), of the form <l>3i,32(k, q) = 
-V2SJ{2 + a) sin 2(9 l>3i,32(k, q), where 



53l(k, q) = 7k(^Xk + Vk)(^q^k-q+Q + ^q^k-q+q) 
+ 7q(^q + ^q)(^k^k-q+Q + ^k^k-q+p) 
+ 7k-q+Q(^k-q+Q + ^k-q+Q)(^k^q + ^k^q) , 

532(k, q) = 7k(^ik + ^k)(^q^k-q+Q + %^k-q+Q) (A23) 
+ 7q(^q + ^q)(^k^k-q+Q + ^k^k-q+p) 
+ 7k-q+Q(^k-q+Q + ^k-q+Q)(^k% + ^k^q) • 



The corresponding self-energy corrections therefore read 

S31 (k, u) = 2 J(2 + a) cos^ sin^ ( A24) 

|^3l(k,q)P 



xE- 

^ a; - CJq - CJk-q+Q 

S32 (k, cj) = -2 J(2 + a) cos^ sin^ 

|<i>32(k,q)P 



zO+ 



(A25) 



E 



a; + cjq + cjQ_k-q - ^0+ 



where 0+ is a small positive number. In taking these 
integrals, 0^ was taken as 10~^. 

Then the on-shell correction to the spectrum is. 



4'^ = ^(l + f)-n^2^E 



|^3i(k,q) 



^k - ^q — ^k-q-hQ 



iO+ 



|^32(k,q) 



^k + ^q + ^Q-k-q ' 



iO+ 



(A26) 



so that combining all the contributions together we ob- 
tain the renormalized spin-wave dispersion 



^k = ^k + fe 



(1) 



J2) 



k -fe^+feL , 



J3) 



(A27) 



wher e £k is the linear spin-wave theory energy, given in 



Eq. (|A14|), 5e\ is the correction due to Hartree-Fock 



decoupling of the quartic terms from Eq. (A19), 5 e\^ i s 
the correction due to angle renormalization in Eq. (A22), 



and dsjj^ is the correction due to cubic vertices from 



Eq. (A26). 
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